#ifndef	INC_MATRIX_H
#define	INC_MATRIX_H

#include <ostream>
#include "BaseTypes.h"
#include "Vector.h"
/*********************************************/
namespace CMatrix
{	
template<typename T>
class CMatrix
{
private:
		
	T		m[4][4];
	
public:

	const T operator[]( const sint32& ) const;			// read only operator
	T& operator[]( const sint32& );					// access operator
	void operator=( const CMatrix<T>& );				// assignment
	CMatrix<T> operator+=( const CMatrix<T>& );			// addition
	void operator*( const T& );					// Multiplication by scalar
	void initModelViewMatrix( CVector::CVector<T>& x, 
				  CVector::CVector<T>& y,
	  			  CVector::CVector<T>& z );		// Basic modelview ( no inversion )
	CMatrix<T> operator!();						// Invert
	
	CMatrix();
	CMatrix( const CMatrix<T>& );
	~CMatrix();		
};
/*********************************************/
template<typename T>
CMatrix<T> CMatrix<T>::operator+=( const CMatrix<T>& mat )
{
	for ( int i = 0; i<4; i++ ) {
		for ( int j = 0; j<4; j++ ) {
			m[i][j] += mat.m[i][j];
		}
	}
}
/*********************************************/
template<typename T>
T& CMatrix<T>::operator[]( const sint32& index )
{
	return( m[index/10][index%10] );
}
/*********************************************/
template<typename T>
const T CMatrix<T>::operator[]( const sint32& index ) const
{
	return( m[index/10][index%10] );
}
/*********************************************/
template<typename T>
void CMatrix<T>::operator=( const CMatrix<T>& mat )
{
	for ( sint32 i = 0; i<4; i++ ) {
		for ( sint32 j = 0; j<4; j++ ) {
			m[i][j] = mat.m[i][j];
		}
	}
}
/*********************************************/
template<typename T>
CMatrix<T>::CMatrix()
{
	for ( sint32 i = 0; i<4; i++ ) {
		for ( sint32 j = 0; j<4; j++ ) {
			i==j ? m[i][j] = 1 : m[i][j] = 0;	
		}
	}
}
/*********************************************/
template<typename T>
CMatrix<T>::CMatrix( const CMatrix<T>& mat )
{
	for ( sint32 i = 0; i<4; i++ ) {
		for ( sint32 j = 0; j<4; j++ ) {
			m[i][j] = mat.m[i][j];
		}
	}
}
/*********************************************/
template<typename T>
CMatrix<T>::~CMatrix()
{
}
/*********************************************/
template<typename T>
void CMatrix<T>::initModelViewMatrix(  	CVector::CVector<T>& x, 
					CVector::CVector<T>& y,
		 			CVector::CVector<T>& z )
{
	m[0][0] = x[0];
	m[1][0] = x[1];
	m[2][0] = x[2];
	m[3][0] = 0;
	
	m[0][1] = y[0];
	m[1][1] = y[1];
	m[2][1] = y[2];
	m[3][1] = 0;
	
	m[0][2] = z[0];
	m[1][2] = z[1];
	m[2][2] = z[2];
	m[3][2] = 0;

	m[0][3] = 0;
	m[1][3] = 0;
	m[2][3] = 0;
	m[3][3] = 1;
}
/*********************************************/
template<typename T>
CMatrix<T> CMatrix<T>::operator!()
{
	CMatrix<T>		temp, inverse;
	sint32			r,c,ri,ci;
	
	// Calculate individual cofactors
	for ( r = 0; r<4; r++ ) {
		for ( c = 0; c<4; c++ ) {
			for ( ri = 0; ri<4; ri++ ) {
				for ( ci = 0; ci<4; ci++ ) {
					temp.m[ri][ci] = m[ri][ci];	
				}
			}
			
			for ( ri = 0; ri<4; ri++ ) {
				temp.m[ri][c] = 0;
			}
			
			for ( ci = 0; ci<4; ci++ ) {
				temp.m[r][ci] = 0;
			}
			
			temp.m[r][c] = 1;
			inverse.m[c][r] = Det(temp);
		}
	}
	
	return( inverse );
}
/*********************************************/












/*********************************************/
//##### Non-member operations
/*********************************************/
template<typename T>
const T Det( const CMatrix<T>& m )
{
	return( 	  m[0]* ( m[11]*( m[22]*m[33] - m[32]*m[23] ) 
			- m[21]*( m[12]*m[33] - m[32]*m[13] )
			+ m[31]*( m[12]*m[23] - m[22]*m[13] ))
			- m[10]*( m[1]*( m[22]*m[33] - m[32]*m[23] )
			- m[21]*( m[2]*m[33] - m[32]*m[3] )
			+ m[31]*( m[2]*m[23] - m[22]*m[3] ))
			+ m[20]*( m[1]*( m[12]*m[33] - m[32]*m[13] )
			- m[11]*( m[2]*m[33] - m[32]*m[3] )
			+ m[31]*( m[2]*m[13] - m[12]*m[3] ))
			- m[30]*( m[1]*( m[12]*m[23] - m[22]*m[13] )
			- m[11]*( m[2]*m[23] - m[22]*m[3] )
			+ m[21]*( m[2]*m[13] - m[12]*m[3] )));
}
/*********************************************/
template<typename T>
CMatrix<T> operator*( const CMatrix<T>& m1, const CMatrix<T>& m2 )
{
	CMatrix<T>		mult;
	sint32			r10;
	
	for ( sint32 r = 0; r<4; r++ ) {
		for ( sint32 c = 0; c<4; c++ ) {
			r10 = r*10;
			mult[r10+c] = 	(m1[r10]*m2[c]) +
					(m1[r10+1]*m2[10+c]) +
					(m1[r10+2]*m2[20+c]) +
					(m1[r10+3]*m2[30+c]);
		}
	}
	
	return( mult );
}
/*********************************************/
template<typename T>
CVector::CVector<T> operator*( const CMatrix<T>& m, const CVector::CVector<T>& v )
{
	CVector::CVector<T>		result;
	sint32				row=0;
	
	for ( sint32 r = 0; r<4; r++ ) {
		row = r*10;
		result[r] = 	m[row]*v[0] + m[row+1]*v[1] +
				m[row+2]*v[2] + m[row+3]*v[3];
	}
	
	return( result );
}
/*********************************************/
template<typename T>
Point3D operator*( const CMatrix<T>& m, const Point3D& v )
{
	Point3D				result;
	fl32				x, y, z, rx, ry, rz;

	v.getCoordinates( x,y,z );
	rx = m[0]*x + m[1]*y + m[2]*z + m[3];
	ry = m[10]*x + m[11]*y + m[12]*z + m[13];
	rz = m[20]*x + m[21]*y + m[22]*z + m[23];
	result.setCoordinates( rx,ry,rz );

	return( result );
}
/*********************************************/
template<typename T>
std::ostream& operator<<( std::ostream& ostr, const CMatrix<T>& m )
{
	using std::endl;
	
	ostr << "[ " << m[0] << " \t" << m[1] << " \t"  << m[2] << " \t"  << m[3] << " ]" << endl;
	ostr << "[ " << m[10] << " \t" << m[11] << " \t"  << m[12] << " \t"  << m[13] << " ]" << endl;
	ostr << "[ " << m[20] << " \t" << m[21] << " \t"  << m[22] << " \t"  << m[23] << " ]" << endl;
	ostr << "[ " << m[30] << " \t" << m[31] << " \t"  << m[32] << " \t"  << m[33] << " ]" << endl;
	
	return( ostr );
}
/*********************************************/
};		// namespace CMatrix
#endif	// INC_MATRIX_H